EPJ manuscript No. 

(will be inserted by the editor) 



^ Thermal critical behavior and universality aspects of the 

o 



CD 



O 

CD 

B 



o 



> 

oo 

o 

O 

(N 
O 
00 
O 



X 



three-dimensional random-field Ising model 



A. Malakis'' and N.G. Fytas^ 

Department of Physics, Section of Solid State Physics, University of Athens, Panepistimiopolis, GR 15784 Zografos, Athens, 
Greece 

Received: date / Revised version: date 

Abstract. The three-dimensional bimodal random-field Ising model is investigated using the N-fold version 
of the Wang-Landau algorithm. The essential energy subspaces are determined by the recently developed 
critical minimum energy subspace technique, and two implementations of this scheme are utilized. The 
random fields are obtained from a bimodal discrete (±/i) distribution, and we study the model for various 
values of the disorder strength A, A — 0.5, 1, 1.5 and 2, on cubic lattices with linear sizes L = 4 — 24. 
We extract information for the probability distributions of the specific heat peaks over samples of random 
fields. This permits us to obtain the phase diagram and present the finite-size behavior of the specific 
heat. The question of saturation of the specific heat is re-examined and it is shown that the open problem 
of universality for the random-field Ising model is strongly influenced by the lack of self-averaging of the 
model. This property appears to be substantially depended on the disorder strength. 

PACS. PACS. 05.50-fq Lattice theory and statistics (Ising, Potts, etc.) - 64.60.Fr Equilibrium properties 
near critical points, critical exponents - 75.10.Nr Spin-glass and other random models 

1 Introduction Hamiltonian of the system is: 

The random-field Ising model (RFIM) [1] is one of the - J S,Sj - ^ /i.S',, (1) 

<ij> i 

most studied glassy magnetic models |2I3I4I5) , mainly be- 
where the Si = ±1 are Ising spins, J is the interaction 

cause of its interest as a simple frustrated system. The 

energy between nearest neighbors, which we take to be 
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e-mail: nfytas@phys.uoa.gr positive so that the model is ferromagnetic and hi are the 
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random fields (RF's). In this paper the values hi are taken transition in 3D exists, provided that the randomness is 

from a bimodal distribution of the form: sufficiently small {Ac ~ 2.3). 

r,/, \ If/L A\ I ^ s^u I A\ However, agreement over several fundamental issues 

is missing and the characterization of the phase transi- 

with A the disorder strength, also called randomness of tion is stiU unclear yL4j. Despite the fact that most stud- 

the system. Various different RF probability distributions jgs suggest a second-order transition [1411511611711511^ . 

have been studied in the past |6|7|8|9|, such as the Gaus- there are also indications of first-order or hybrid-order 

sian distribution, the wide bimodal distribution (with a transition [HID]. Note also that the mean field theory [H] 

Gaussian width), and the bimodal distribution considered differentiates between a binary and a continuous random- 

also here. ness distribution, predicting for the former a tricritical 

In spite of many years of study, the critical behavior point at which the transition becomes of the first-order, 

of the three-dimensional (3D) RFIM has been a matter of at high fields. However, it is now generally accepted that 

several debates and is still controversial. One of the early a new fixed point controls the behavior of RF ferromag- 

disagreements was the question of whether the model un- nets [22123] . The significance of this for the RFIM (in d > 

dergoes a phase transition from a high temperature para- 2) is that this new zero temperature random fixed point 

magnetic phase to a low temperature ferromagnetic one, controls the whole critical line (TdA)) and that the RF's 

for some range of the randomness A. The work of Parisi are always relevant. For disordered systems with weak 

and Sourlas 10 introduced the notion of dimensional re- randomness which couples to the local energy (such as 

duction, indicating that the critical behavior of the RFIM random-site impurity or random-bond models) the crossover 

in d dimensions, at sufficiently low randomness, should be to a new random fixed point, depends on the Harris cri- 

identical to that of the well-known normal Ising model in terion [23124] . According to this, the disorder is relevant 

d — 2 dimensions. This in turn indicated that the model if the correlation length exponent of the pure model {v — 

should not exhibit a phase transition in 3D or fewer. How- impure) satisfies the condition dv < 2 and this condition 

ever, a different argument based on the droplet theory of may be stated, with the help of the hyper-scaling relation 

domain wall energies in the ferromagnetic state [11] , sug- {a — 2 — dv), as a > 0. Since the specific heat exponent 

gested that a phase transition should exist in 3D for finite of the 3D Ising model is positive, weak disorder should be 

temperature and randomness. The whole puzzle has been expected to be relevant. In the case of the RFIM the type 

largely cleared out by Imbrie [12] and Bricmont and Kupi- of disorder is much more severe, since the randomness cou- 

ainen [13| . who showed the existence of an ordered phase, pies to the local order-parameter and the crossover renor- 

Their arguments strongly supported the view that a phase malization group eigenvalue is always positive [23| . The 
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inequality v > 2/d, derived by Chayes et al. [22] for the A relevant active and enigmatic issue concerns the be- 
correlation length exponent of a generic disordered system havior of the specific heat (see Ref. [29] and references 
{v — Vrmidom) would imply, using again hyper-scaling, a therein). The specific heat of the RFIM can be experimen- 
negative specific heat exponent {a < 0). However, it is tally measured and is of considerable theoretical interest, 
believed that hyper-scaling is violated in the RFIM and There is a strong disagreement in literature about the pos- 
the specific heat exponent a is related to by a modified sible divergence or saturation of the specific heat. In stud- 
hyper-scaling law 2 ~ a — (d — O)^, where the exponent ies supporting the scenario of saturation there is a discrep- 
9 characterizes the scaling of the stiffness of the ordered ancy in the reported negative values of the critical expo- 
phase at the critical point. Thus, the specific heat expo- nent a. Some of these later studies find strongly negative 
nent of the RFIM is not restricted, by the above theoret- values, ranging from a — —1.5 [3T] to a = —0.5 |6I14I32| . 
ical considerations, to be negative |22j . In particular, Hartmann and Young ^ recently found by 

a ground state technique the value a = —0.63 ± 0.07, 



whereas Middleton and Fisher [33|, using the same tech- 
nique, estimated in marked disagreement a — —0.01 ± 
0.09. 



A general sketch of the phase diagram of the RFIM is 
given in several papers [618125] and will be also presented 
here in Sect. 13.31 At low temperatures and moderate values 

of randomness, the system is assumed to be in an ordered From the experimental point of view, a true realiza- 

ferromagnetic phase, whereas in the opposite regime the tion of the RFIM is hardly conceived. However, it has 

system is paramagnetic. From the notions of the perturba- been shown that dilute antiferromagnets in uniform ex- 

tive renormalization group (PRG) it is expected that the ternal field (DAFF) represent physical realizations of the 

RF is the relevant perturbation at the pure {A = 0) fixed RFIM [34] and a number of experiments investigated the 

point, and that the RF fixed point is at T = 0. However, it phase transitions of such 3D systems [3^. These experi- 

is known that PRG fails for the RFIM and that a theoret- ments have proven to be very difficult and their interpreta- 

ical justification of universality for this and also other dis- tion doubtful due to the extremely slow, glassy dynamics 

ordered systems is lacking |2I3I9I10| . Questions concerning of the system. Experiments on DAFF, provided evidence 

the general characterization of the phase transition, the of a second-order phase transition and a logarithmic singu- 

existence of an intermediate glassy phase [25126127] , the larity for the specific heat [3^ . Note that recently. Barber 

behavior of the renormalization group flow in the middle and Belanger [39] in their Monte Carlo study of a DAFF 

of the phase diagram _25J, and the dependence of the criti- model reported also that their specific heat curve closely 

cal exponents on the randomness distribution and disorder mimics a logarithmic peak. Their results are based on a 

strength are still open [719128] . large lattice (L = 128) but instead of sample averaging 
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they have observed the behavior of only a few RF real- 
izations. On the other hand, there is also experimental 
evidence |40| supporting the opposite view of a cusp-like 
singularity of the specific heat, in agreement with a sat- 
urating specific heat (a < 0) as found in the studies of 
Refs. [Hn]. 

It has been pointed out that a strongly negative value 
of a causes serious difficulties with respect to the Rush- 
brooke relation: a + 2/3 + 7 > 2 |6|14|32|33j . Therefore, 
there have been several attempts |6|15|16ll7|41j in order 
to find a consistent set of scaling relations to describe the 
critical behavior of the RFIM. Among the several scaling 
scenarios proposed |9|26|27|28|32|33j . the single second- 
order critical point behavior characterized by three scaling 
exponents [33j seems to be consistent with a close to zero 
estimate for the specific heat exponent. Thus, the above 
described conflicting situation in literature concerning the 
divergence or saturation of the specific heat is one of the 
open important issues, whose implications on the critical 
behavior of the model are not understood. 

The first step towards its resolution was taken recently 
by the present authors [29l30j . where an extended nu- 
merical investigation of the 3D bimodal (A = 2) RFIM 
revealed the importance of the property of lack of self- 
averaging of the specific heat of the model, as well as the 
possibility of large-L crossover phenomena in the scaling 
behavior of the specific heat of the model. Here, we extend 
our analysis for the values Z\ = 0.5, 1 and 2 of the disor- 
der strength below the critical value Ac in order to obtain 
a more comprehensive picture. To this end, we implement 
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recently developed efficient Monte Carlo methods that di- 
rectly calculate the density of states (DOS) of a classi- 
cal statistical model. A brief overview of the numerical 
techniques used in the past for the RFIM are presented 
in the next Section, together with the necessary details 
of the methods used in our approach. The utilization of 
our recently proposed critical minimum energy subspace 
(CrMES) scheme [12] to the RFIM is also explained. In 
Sect.|3]the new numerical results for the cases A — 0.5, 1 
and Z\ = 1.5 are given and the phase diagram of the model 
is reproduced. The universality aspects of the model are 
discussed and found to support the scenario of violation of 
universality. Our conclusions are summarized in Sect.|4l 

2 Numerical techniques 

There are two distinct kinds of numerical approaches for 
the RFIM. In the first approach, traditional Monte Carlo 
methods are used to simulate the properties of the system 
at finite temperatures |14I18I31I32I39I43I44I45| . The sec- 
ond approach is grounded on the well-known belief that 
the critical behavior of the model is governed by the non 
trivial RF fixed point at T = [33146] . In this case, 
graph theoretical algorithms |7l47j are used to calculate 
the ground states of the system for a sample of RF's at 
different disorder strength. Using this later approach quite 
large lattices have been studied: L < 80 \TM , L < 90 [S], 
L < 96 [6 and finally L < 256 [33]. Yet, in the traditional 
Monte Carlo approach the sizes studied were restricted to 
the size L < 16 |8ll4l31j . L < 20 [32] and finally we may 
refer, as an exception, to the case L — 128 in the study 
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of Ref. [39] for particular RF's as mentioned in the intro- advance. To obtain a good approximation of the locations 
duction. From the T = numerical studies one obtains of the specific heat peaks, the temperature step must be 
an accurate estimate of the critical randomness and from chosen sufficiently small for, otherwise any interpolation 
the finite temperature studies further information for the scheme may miss the correct height of a possible sharp 
phase diagram may be derived. From the finite tempera- peak. In fact, this situation of a possible sharp peak, turns 
ture approach one can also find, by extrapolation, a crude out for a significant number of RF's [29j . 
estimate for the critical randomness [S] (see Sect. [3]). It 
is worth noting that quite recently Wu and Machta, com- 
bining finite and zero temperature studies of the RFIM, 
reported strong correlations of ground states and thermal 
states near the critical line for given realizations of the 
disorder, supporting strongly the T = fixed point sce- 
nario l48l. 



From the above discussion one should wonder whether 
the traditional Metropolis sampling could be trusted to 
provide even a moderate estimation plan, since it requires 
immense computer resources and faces all mentioned prob- 
lems. The cluster flipping algorithm for the RFIM pro- 
posed by Dotsenko, Selke and Talapov [IHj is a straight- 
In traditional Monte Carlo studies of the RFIM |8l31ll8ll41^ ard extension of the Wolff algorithm [50], devised to 
the system is simulated in a restricted range of temper- overcome the slowing down effect and speed up the flip 
atures, appropriate for the location of the pseudocritical dynamics. A more efficient form of this algorithm, the 
temperatures. However, single spin-flip methods, such as limited cluster flip (LCF) algorithm, has been invented 
the Metropolis or the heat bath algorithms, face severe by Newman and Barkema [5] and was used for the study 
slowing down problems of equilibration and temperature of the Gaussian RFIM. Furthermore, these authors have 
averaging since the characteristic times may be exponen- combined the LCF algorithm with the generalized his- 
tially large at low temperatures (T < Tc), as explained togram method of Ferrenberg and Swendsen [51] which 
in Ref. |8J. Moreover, the sample averaging process in- is a re- weighting scheme, using a restricted set of temper- 
troduces new characteristic features and requires further ature measurements. This combination may be hopefully 
computer resources. Indeed, the appropriate pseudocrit- more reliable for the location of the pseudocritical tem- 
ical temperature for the RFIM is a strongly fluctuating peratures. Finally, a new cluster technique, that combines 
quantity [29] . and this property amplifies the computer the replica-exchange method of Swendsen and Wang |52| 
time requirements for its location. Hence, depending on and the two- replica cluster method [53j . was implemented 
the size of the lattice and the disorder strength, it is nec- by Machta, Newman, and Chayes [33] where single real- 
essary to simulate the system for each RF realization in a izations of the disorder strength were studied for sizes up 
quite wide temperature range, which is not even known in to L — 24^^. 



6 A. Malakis and N.G. Fytas: Thermal critical behavior and universality aspects of the three-dimensional... 

Here, we employ a different strategy which utilizes the trial. Here, of course, fj is the value of the WL modi- 
new and popular methods of efficient estimation of spec- fication factor / |59| at the jth iteration, in the process 
tral degeneracies of classical statistical models |44l51l55l56l57(3'8l&DI6l0lDll68f reducing its value to 1 , where the de- 
and the recently developed CrMES technique [42]. This tailed balance condition is satisfied. In all our simulations 
scheme has the merit of locating the pseudocritical tem- the WL modification or control parameter takes the ini- 
peratures by determining the DOS in the proper energy tial value: = e « 2.71828.... When starting a new 
subspace by using simple algorithms in a unified imple- iteration the control parameter is changed according to 
mentation. Moreover, it avoids all the above problems, fj+i = \/7j,j = 1,2,..., 20 [59]. For the histogram flat- 
speeding up the simulations. Specifically, we use the multi- ness criterion, we use a flatness level 0.05, as in previous 
range Wang-Landau (WL) algorithm [55], and its N-fold studies |42I55] . 
version as presented by Schulz et al. |60| . The accuracy of 

this scheme was discussed in Ref. [2^, where more details 2.2 The N-fold version of the Wang-Landau algorithm 



than those given below for the appropriate implementa- 
tions can be found. 



For the bimodal RFIM is convenient to use an index n 
to characterize directly the corresponding energy changes 

produced by the spin-flip process. The number of differ- 

2.1 The Wang-Landau algorithm . , r at r i i i kt ^ i 

° ° ent classes tor the N-fold version n — 1,...,A/ depends 



on the value of the disorder strength A. For example, 
consider the case A ^ 1. The possible energy changes 
are AEn — ±14, ±10, ±6 and ±2, and using an in- 
dex ri = 1, 2, 8 corresponding to 8 classes we can write 
AEn = -14+(n-l)-4. Note that, the RF value at the site 
in which the spin-flip is going to take place is also affecting 
the energy change. Denoting the populations of spins by 
Nn, J2n -^n — N, whcrc N is the number of sites: N — L^, 
the selection probability of a class, P„, will be propor- 
^ ^ min{l, G{E)/G{E + AE)}, {E + AE) e I ^j^j^^^j ^-^lis number multiplied by the corresponding ac- 

0, {E + AE) ^ I ceptance ratio. For the application of the algorithm in 

the random walk is not allowed to move outside of the en- multi-range approach we follow the description of Schulz 
ergy range, and we always increment the histogram H{E) — » et al. 60J for the N-fold version of the WL method. If the 
H{E) + 1 and the DOS G{E) — *■ G{E) * fj after a spin-flip system is in a spin state with energy E I and after the 



For the application of the WL algorithm in a multi-range 
approach we follow the description of Schulz et al. [60] , i.e. , 
whenever the energy range is restricted we use the updat- 
ing scheme 2 in that paper. Consider the restriction of the 
random walk in a particular energy range / = [Ei,E2] 
and assume that the random walk is at the border of the 
range /. Then, the next spin-flip attempt is determined 
by the modifled Metropolis acceptance ratio: 
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spin- flip is in a state with energy E' = E + AEn, the se- where {E^,E+) is the minimum dominant subspace sat- 
lection of the class for the next spin-flip is obtained from isfying the following accuracy criterion: 
the following [55] : 



Cl{E-,E+) ^ 



r (F F ^ ^ - 

min{l, G(i?)/G(£")}, E' E I in this paper we have used the accuracy criterion r = 

0, E' ^ I lO"'*, which is extremely demanding compared to the rel- 

(4) 

ative errors produced in the specific heat, say by the WL 

The average life time of a state, reflecting the number of 

method. It is also a very strict criterion for the present 

attempts we expect the system to remain in its current 

model, in view of the existing very large sample-to-sample 

state before moving to the new state, is At — W/Z [BD] 

fluctuations of the specific heat. A practical algorithmic 

where: 

approach for specifying the CrMES is fully described in 
^J-' " ^ Ref. [H]. We may satisfy the specific heat accuracy crite- 

The rest of details for the algorithms can be found in the ria defined in equation ^ for any particular realization 
original papers [55160] . of the RF, by restricting the WL random walk in the cor- 

responding critical energy subspace. 

2.3 Implementation of the CrMES technique ^his restriction greatly facilitates our simulations with- 

out introducing additional errors. Since we don't know in 
The CrMES scheme [42] uses only a small part {E-,E+) advance the CrMES for a specific rcahzation of the RF 
of the energy space {E^in.Emax) to determine the spe- ^e have two alternatives. The first option is to use an efii- 
cific heat peaks. If E is the value of energy producing the cient prognostic method of identifying the CrMES for any 
maximum term in the partition function at the tempera- particular realization of the RF by using the first stages 
ture of interest (say the pseudocritical temperature), the of the WL method. For instance, one may try to estimate 
sums are restricted as follows: the CrMES from the first 12 iterations in the process of re- 

ducing the WL modification factor /. To implement safely 
this option, one should be careful to use for the rest of 

WL iterations a much wider energy range than that pre- 

(6) . . 

dieted in the first 12 iterations. The second option is to 
'guess' (by using some convenient extrapolation method) 

and 

~ a broad energy subspace that will cover the overlap of the 

<^{E) = [S{E)-I3E]-\S{E) - fiEy, ^ = ^ [^(^)]. CrMES for all RF's in the sample. Implementing the sec- 

(7) ond method is straightforward and has the advantage that 
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the approximation for the specific heat curve of a particu- aged specific heat, without raising the question of whether 

lar RF realization is accurate in a wide temperature range this averaged curve is the proper statistical representative 

including the pseudocritical temperature corresponding to of the system. The peak of this averaged curve was then 

the particular RF. This option (the second) was used for analyzed by using finite-size scaling relations. On the other 

the cases A = 0.5 and 1, whereas for A — 1.5 (and 2 29J) hand, the work of Barber and Belanger [39j, in which the 

both options were used, each for 50% of the simulations behavior was observed for particular RF realizations, is 

(for a more detailed discussion see Ref. :29J). a different route but it would be hard to accept this as 

an adequate alternative. A particular RF is not generally 

_ _ , , , . representative of the behavior of a large sample of RF's. 

3 Results and analysis 

Indeed, in previous papers [611413113^ the following 
3.1 Definitions and the property of self-averaging average has been considered for the specific heat: 

M 



For a disordered system we have to perform two distinct l- j"." ^ 

m— 1 

and the finite-size scaling behavior of the peak of this 

usual thermal average has to be carried out and secondly 



kinds of averaging. Firstly, for each RF realization the 

averaged curve has been studied, by assuming that the 

we have to average over the distribution of the random . rr^i i i i i- 

maximum [CJ*„ — ma.XT{[C\avs and the correspondmg 

parameters. In effect, this means that we must gener- , . . , , , , ,. , 

pseudocritical temperature obey the scafmg faws: 

ate and study quite large samples of RF's. Following the 

[C]L =P + cL"/'^; Tl{[C]av)=T^ + bL-^/'', (10) 

methods outlined in Sect. 12.31 the thermal average for the 



specific heat is given by the approximations ©-(IH]). Us- where a and ly are considered to be the specific heat and 
ing large samples of RF's we can estimate the relevant correlation length critical exponents, respectively. Note 
probability distributions of the location of the specific that, these averaged curves are very sensitive to the prop- 
heat peaks. Let C„i(T) denote the specific heat of a par- erty of lack of self-averaging (see the discussion below) 
ticular realization m in a sample of M realizations of due to the fact that the corresponding thermodynamic 
the RF (m = 1,2,..., M). The pseudocritical temperature quantities are characterized by broad distributions in the 
Tl{Cm{T)) depends on the realization of the RF and for thermodynamic limit [55]. 

large values of the randomness A, is a strongly fluctuat- It is clear that when studying random systems the only 
ing quantity [29 . Let us also denote the locations of the meaningful objects for investigating the finite-size scaling 
specific heat peaks by (Cj^,T£ ^). It seems that, in all behavior are the distributions of various properties in en- 
previous studies |6ll4l31j . the averaging process over an sembles of several realizations of the randomness. Hence, 
ensemble of RF's was carried out on the curve of the aver- it is important to be able to ascertain to what extent are 
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the results obtained from an ensemble of random realiza- 
tions representative of the general class to which the sys- 
tem belongs. The answer hinges on the important issue 
of self-averaging. In a self-averaging system, a single very 
large system suffices to represent the ensemble; without 
self-averaging, a measurement performed in a single sam- 
ple, no matter how large, does not give a meaningful result 
and must be repeated on many samples. In a Monte Carlo 
study of a self-averaging disordered system the number 
of samples needed to obtain the average [Q] (e.g., Q can 
be the energy, magnetization, specific heat, or susceptibil- 
ity) to a given relative accuracy decreases with increas- 
ing L. On the other hand, in a non self-averaging system, 
the number of samples that must be simulated rises very 
strongly with L. If a quantity is not self-averaging, then 
we talk about lack of self-averaging and as explained the 
process of increasing L does not improve the statistics. 
In other words, the sample-to-sample fluctuations remain 
large. The problem of self-averaging in the RFIM has been 
a matter of investigation over the last years [29130163] . A 
common measure characterizing the self-averaging prop- 
erty of a system based on the theory of finite-size scaling 
has been discussed by Binder [Blj and has been used for 
the study of some random systems |65I66] . This measure 
inspects the behavior of a normalized square width quan- 
tity, defined as: 

Vq 



Rq 



(11) 



where Vq = [Q^] — [Q]^ is the samplc-to-sample variance 
of the average [Q]. Here, Q is used in respect of the specific 
heat Cj^. According to the literature [64|65|66] when the 
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Fig. 1. L-dependence of the ratio -Rc*, defined in equation (lllfl 
for A = 0.5, 1 and 1.5. The inset illustrates the variation of 
Rc^ a-s a function of the disorder strength A, including the 
value for the case A = 2 [29] . 

ratio Rq tends to a constant, the system is said to be 
non self-averaging and the corresponding distribution (say 
P{Q)) does not become sharp in the thermodynamic limit. 

In Rcf. 29J it has been shown that the specific heat of 
the bimodal RFIM for the case Zi = 2 is characterized by 
the property of lack of self-averaging (see inset of figure 4 
in Ref. [IH])- In analogy with the case Z\ = 2 of Ref. ^5] . 
we construct the ratio Rc*^ for the cases A — 0.5, 1 and 
1.5 and plot the results in figure [TJ The data presented 
for the cases A — 0.5 and 1 are taken from 1000 samples 
of RF's for i < 10 and 400 for L = 14 - 24, while for 
the case A — 1.5 the samples of RF's used were: 1000 for 
i < 12 and 300 for i = 12 - 24. 

From figure [T] we observe that for small values of the 
randomness A the ratio Rc:^ has values close to zero. In 
particular, for A — 0.5 the ratio Rq* shows a rather faint 
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dependence on L and is practically very close to zero. For 
this case, and possibly for even weaker RFs, it appears 
that Rci^ ~^ and the property of self-averaging may be 
well obeyed. However, we know from our previous study 
for Z\ = 2 that the above property is strongly violated, at 
the same lattice sizes, and that Rc^^ 0-3 29J. Thus, for 
strong disorder the behavior appears very different and 
the development of the increasing influence of the ran- 
domness A on the ratio Rc^ can be seen in figure [T] and 
in particular from the corresponding inset. For the lattice 
sizes studied here (L = 4 — 24) the ratio Rc^ for the cases 
A > 1 seems to increase with the lattice size and the es- 
timated non-zero limiting values for A = 1 and 1.5 are 
Rc^ 0.016 and Rc^ ~^ 0.12, respectively. However, 
the behavior of the specific heat is notoriously difficult 
even in simple pure models |64| , and for the present model 
the already existing conflicting situation is an additional 
warning against drawing definite conclusions at these lat- 
tice sizes. It is quite possible that we are not yet in the 
regime of large enough L, where simple scaling laws would 
be expected to hold. For the case Z\ = 2 we have already 
observed [53] that the system appears to crossover and 
change behavior for sizes L > 32 and we have suggested 
in that paper that the finite-size study should be extended 
to at least L = 80 in order to have a more convincing pic- 
ture. For the strengths studied here, A — 0.5, 1 and 1.5, 
we suspect that even much larger sizes would be needed in 
order to draw definite conclusions for the true asymptotic 
behavior. 





3,0 p 




2,7 - 


> 

CD 


2,4 - 


O 




"> 

CD 


2,1 - 


E 




o 


1,8 - 




1,5 - 




1,2 - 




Fig. 2. Finite-size behavior of [C^]av and [C']av for A = 0.5. 
The solid and dotted lines correspond to logarithmic fits for 
[C'^]av and [C]*„, respectively. 



There are several cases in the literature where the char- 
acterization of a phase transition demands very large lin- 
ear sizes and the picture obtained from moderate sizes is 
completely misleading. A characteristic example is the 5- 
state 2D Potts model, for which Landau [S7] suggested 
that the expected first-order behavior would not be clar- 
ified from finite-size data up to sizes L = 2000. In a dif- 
ferent inquiry Hilfer et al. |68| estimated that the asymp- 
totic behavior of the tail regime of the universal order- 
parameter distribution for the 2D Ising model would re- 
quire sizes of the order of L > 10^. Thus, having to deal 
with the controversial 3D RFIM, for which even the ex- 
istence a tricritical point at high fields is not yet clari- 
fied [69j . we prefer to regard the observed in figure [T] strong 
violation of the self-averaging property as a rather tenta- 
tive conclusion which has to be further verified by studying 
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larger systems and more physical properties (such as the 
magnetic susceptibility [30]). 

Since all past finite temperature studies were attempted 
on small and moderate sizes (L < 20), it is valuable to ex- 
amine the implications of the strong violation of the self- 
averaging property at these moderate sizes. In our opinion 
the inconsistent estimations in the literature have, at least 
partly, their origin on such an unconventional behavior of 
the RFIM. In order, to observe better these implications 
we proceed to study, in addition to the above scaling laws, 
the sample averages of the individual specific heat maxima 
and pseudocritical temperatures defined by: 
1 



\C1 



M 
1 



L,m\av 



J^J2Km = T, + bL-'/r (12) 



The mean values defined above characterize the corre- 
sponding probability distributions and consist a different 
kind of representative of the samples of RF's. To quan- 
tify the sample-to-sample variations of the specific heat 
peaks we use the standard deviation of C*j over a sam- 
ple of m — 1,2,...,M RF realizations, Vc^- This is the 
parameter of equation (jlip and figure [T] and will be also 
illustrated in the following figures as error bars. However, 
it should not be in any case confused with the statistical 
errors resulting from the thermal average approximations 
of equations ^ and ([7]) . 

3.2 Scaling behavior of the specific heat 

Let us start by presenting in figure [5] the finite-size be- 
havior of the peaks of the sample average [C^]av and that 




Fig. 3. Finite-size behavior of [C^]av and [C]av for A — 1. The 
solid and dotted lines correspond to power law fits for [C^lau 
and [C] Jt,, respectively. Both quantities saturate, however with 
different exponents. 

of the averaged curve [C]*^, for the case A = 0.5. The 
number of RF realizations is M — 1000 for L < 10 and 
M ^ 400 for L = 12 - 24. The difference between the be- 
havior of the peaks of the averaged curve [C]*^ and that of 
the sample average [C^]av does not emerge for small val- 
ues of the lattice size L. Only for L > 16 is the sample-to- 
sample fluctuation considerable and seems to differentiate, 
although mildly, between the two averages. Noteworthy 
that, in this case the standard deviation of the sample- 
to-sample fluctuations is signiflcantly smaller than the av- 
erage [C;^]av- Vc;^ < [C]*„ < [C,*Jq.u. Based on the data 
L = 6 — 24, no sign of saturation for both [Cj^Jan and [C]*„ 
is observed, and one would be tempted to predict a mildly 
diverging behavior. In fact the best flts, corresponding to 
the smallest value of P^r degree of freedom, predict a 
logarithmic behavior of the form [C*Jau = 0.844(3) • InL 
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and [C]*„ = 0.834(4) • InL, respectively. Note that, the fit 
for [Cm]at) has a smaller value of per degree of freedom 
than that of the fit for [C]*„. Attempting a power law for 
[G^]av we find a diverging behavior with a much larger 
value of per degree of freedom. 

Next, we consider the intermediate case where the ran- 
domness takes the value A = 1. Figure [3] shows again 
the finite-size behavior of the peaks of the sample average 
[C^]av and that of the averaged curve [C]*^. In this case 
we apply power law fits of the form: [C]*„ — P + cL^ 
and [C*Jau = p + cL^ , as the one proposed in equa- 
tions (|10p and using the same number of RF real- 
izations as in the case A — 0.5. These fits yield saturation 
laws for both cases [C*Jat) and [C]*^, predicting however 
different saturation values 4.15(65) and 2.88(13), respec- 
tively. Specifically, the relevant fits of comparable quan- 
tity, give: {p,c,w) = (4.15(65), -4.68(44), -0.31(9)) and 
(p,c,w) = (2.88(13), -3.56(64), -0.51(56)). From these, 
we could even speculate that the saturation of both quan- 
tities takes place with different exponents: w = —0.31 and 
w = —0.51. This value of w corresponding to the peaks of 
the averaged curve [C]*^^ compares very well to the value 
-0.5 of Ref. 31J for the ratio h/T = 0.25 of their study [h 
is used for the disorder strength in Ref. |31|). Using our 
approximate phase diagram (see below figure [5]), we find 
that their case closely corresponds to the case A — 1 stud- 
ied here. However, the values of the exponents estimated 
from the fits can not be taken seriously, since this analysis 
does not account for the systematic problem that at least 
part of the data are not yet in the regime of large enough 




Fig. 4. The same as in figure [31 but for A — 1.5. The clear 
and early saturation of [C]^^ is similar to that of the case 
A = 2 [2S]. The diverging power law behavior of [C^]aD shown 
by the solid line gives an exponent w — 0.44(7). 

L, where finite-size scaling without corrections holds. Note 
that, the standard deviation of the sample-to-sample fluc- 
tuations in this case is also smaller than [Cj^Jai,, that is 
Vc^ < [C]*„ < [C^]av, but not as small as in the case 
A = 0.5. 

Finally we treat the case A ^ 1.5. Figure |4] illustrates 
the finite-size behavior of the peaks of the sample average 
[Cmlau ^^'^ that of the averaged curve [C]*^,. The data 
presented here are taken from samples of M ~ 1000 RF's 
for L < 12 and M = 300 for L = 14 - 24. The satura- 
tion of the peaks for the averaged curve is quite obvious 
and is attained, as in the case A = 2 [29], already in 
the small L-regime. Furthermore, the behavior of the av- 
erage [C*Jat) looks similar with that of the case A = 2 
for L < 24, and despite the fact that there are no signs 
of saturation of this quantity for these lattice sizes, the 
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Fig. 5. Size dependence of pseudocritical temperatures for var- 
ious values of A, including the case Z\ = of the normal cubic 
Ising model [42], and the case A = 2 [29] . 

possibility of crossing over to a final saturation for larger 
lattice sizes can not be excluded. It is worth noting that, 
as in the case A = 2, the standard deviation of the sample- 
to-sample fluctuations seems to obey the same behavior 
with that of [C;^]av, since Vb;; ~ 2{[C;^]av - [C]^^) and 
[C]*^, ~ 1.48. Our power law fitting attempts predict for 
[C^]av a diverging behavior with p = — 0.21(3), 0.7(2), 
and w = 0.44(7). Meanwhile, the averaged curve [C]*„ 
strongly saturates with p = 1.48(2), c ~ —5.15(1.6) and 
an exponent w — —1.69(24), already from the small L- 
regime. The saturation exponent of the averaged curve 
w — ajv = —1.69(24) should be compared to the value 
w — a/iy = —1.1(4) given in Ref. [3T] for h/T = 0.5, which 
now corresponds approximately to our Z\ = 1.5 case (see 
figure [6]). 

Comparing the behavior of [Cj^Jat, for the cases A = 1 
and A = 1.5, one may discern a conflicting picture, in 



a sense that while w is negative for Z\ = 1 - indicating a 
strong saturation - the same exponent turns out to be pos- 
itive - indicating a rather strong divergence - for A = 1.5. 
We believe though, that this is not a surprise. In fact, 
the blowing up of the property of lack of self-averaging 
in the range A = 1 — 1.5, as illustrated in the inset of 
figure [1] may be behind this behavior. A possible satu- 
ration in the asymptotic limit may occur in both cases 
but this may happen via different complex routes because 
of the unsettled and (Z\, i)-sensitive self-averaging prop- 
erty of the system. On the other hand, the quantity [C]*^, 
is very weakly i-depended in the large i-regime and its 
early saturation to a value (that depends on the disorder 
strength), leaves no room for an accurate estimation of 
its behavior since the statistical errors dominate in the 
large L-regime. This fact, when combined with the pos- 
sible crossover behavior of the system at quite large lin- 
ear sizes, larger than those corresponding to the above 
discussed saturation, lead us to suggest that scaling at- 
tempts on [C]*^, including previous studies, should not 
be trusted. 

3.3 Phase diagram and universality aspects 

In figure we present the size dependence of the pseud- 
ocritical temperatures for all values of randomness stud- 
ied. We have included the case of the normal cubic Ising 
model, for which the numerical data of Ref. [42] have been 
used, and the case A = 2 |5S], using results up to i = 24, 
where our numerical scheme is accurate. The results of 
the power law fittings applied (see equation ([T^)) are pre- 
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Table 1. Critical temperatures and exponents for various val- i.e. Refs. [8125] ). The dotted line shows a power law fit of 
ues of A, including the case Z\ = of the pure Ising model (42j, the form: 

and the case Z\ = 2 [29]. ^ p + qA'' , (13) 

with p = 4.5114(20), q = -0.59(3), and r = 2.29(6). The 
value of p is very close to the value of the critical tem- 
perature of the normal 3D Ising model (Tc « 4.51153...), 
approving to a certain extent our fitting choice. The above 
power law ansantz for the phase diagram Tc — Tc{A) 
has a clear physical motivation, which could be compared 



A 




b 







4.51153 


-3.96(14) 


0.66(6) 


0.5 


4.380(3) 


-4.12(51) 


0.57(24) 


1 


3.949(19) 


1.40(83) 


0.95(28) 


1.5 


3.001(13) 


1.85(48) 


1.86(34) 


2 


1.63(19) 


2.28(8) 


2.55(49) 



with various functions considered in the past [8154] , The 

sented in table [T] From table [T] it is obvious that there is a critical value of the randomness is estimated to be = 

strong dependence of the shift exponent on the value 2.42 ± 0.18, close to the value 2.3 ± 0.2 of Newman and 

of the disorder strength. While, for relatively small values Barkema and the value 2.35 of Ogielski [18]. To get a 

of the disorder strength A the shifting of the pseudocrit- better estimate for Ac the system should be simulated for 

ical temperatures follows that of the normal Ising model, larger values of A close enough to the critical value, 
for larger values of A the exponent D shows an intense 
variation, indicating a possible violation of universality, 

4 Concluding remarks 

in agreement with the results of Sour las 9 . In fact, it 

is known that the only theoretical arguments supporting The numerical route utilized here for the study of the 

the existence of universality classes in random systems are RFIM consisted of the application of the multi-range WL 

based on PRG theory and these arguments have been in- algorithm [59 in its N-fold version [50], implemented within 

tensively called into questioned for the case of the RFIM. the CrMES scheme l42j. We hope that the presented com- 

Equivalent studies of universality violations have been re- bination of algorithms and techniques will be useful in 

ported also in other glassy systems [7D], reenforcing the further numerical studies of this and other similarly chal- 

view that the concept of universality in complex systems lenging problems, such as spin glasses, 

is not fully clarified and that more work needs to be done Qur analysis showed that, in general, the behavior of 

towards this direction. the mean [CZ]av is distinct from that of [C]*av and that 

Based on the data of table [2 we give in figure [5] an this is a result of the lack of self-averaging, a property 

approximation of the phase diagram of the model which that varies strongly with the disorder strength and the 

is comparable with the ones given in the literature (see lattice sizes considered. For sufficiently small values of 
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Fig. 6. The phase diagram based on the data of table [T] and 
equation (|13p . The dotted line separates the ferromagnetic (F) 
from the paramagnetic (P) phase. The estimated critical value 
Ac of the disorder strength, above which no phase transition 
occurs, is Ac = 2.42 ± 0.12. 

randomness Zi, [C*Jau and [C]*„ seem to obey a mildly 
diverging behavior, showing no signs of saturation. Mov- 
ing to the intermediate range (A — 1), both quantities 
saturate, but with different exponents. Yet, this area of 
the phase diagram of the RFIM is not fully understood. 
Theoretical studies based on the replica formalism predict 
the existence of an intermediate glassy-like phase which is 
characterized by a breaking of replica symmetry near the 
transition temperature 5 . However, the interpretation of 
these results is not clear, and understanding may be sim- 
pler under the prism of a strong and complex variation 
of the property of lack of self-averaging. For large values 
of randomness, say A — 1.5, the peak of the averaged 
specific heat curve [C]*^ obeys a strong saturation law 
which is attained already in the small L-regime, in agree- 



ment with our previous findings for an even larger value 
of the disorder strength (case A — 2 [29]). But as men- 
tioned earlier the corresponding scaling attempts would be 
hardly trusted. In the same range of the disorder strength, 
the behavior of the sample mean [Cj^Jau is somewhat sur- 
prising, showing no signs of saturation, and its behavior 
seems to follow the large sample-to-sample fluctuations 
developed by the blowing up of the property of the lack of 
self-averaging. Apparently, the blowing up of the property 
of lack of self- averaging in the case = 1.5 is responsible 
for this behavior, shifting a possible saturation to larger 
values of L. Provided that our previous analysis for the 
case A = 2 [29j recorded a 'final and unexpected' saturat- 
ing behavior of the sample average [Cm] aw for L > 32, it 
will be interesting to observe whether this behavior pre- 
vails for larger lattice sizes, even in cases of small values 
of randomness. 

Turning to the shift behavior of the pseudocritical tem- 
peratures of the model, we found a very strong dependence 
of the shift exponent on the disorder strength, reinforcing 
the scenario of universality violation. The shift for small 
values of A appears to follow the direction of the pure 
case (A = 0) of shifting to Tc{A) from below and this 
is reflected in the negative sign of the parameter b. For 
large values of A, the power law exponent ? shows a very 
strong variation, which may be due to the existence of ad- 
ditional leading and non- leading correction terms f9!. In 
order to support numerically the concept of universality 
for the exponent v one should have accurate data for very 
large lattices, as has been pointed out also in previous 
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